function [mergers,subs,maxs,Omega]=equilibrium(sigmaS,sigmaT,N,T,lambda1,lambda2,q_team,q_sp,nScore,Svalues,mergers_data,nSim)

%Values
L_sp_end=1;
L_team_end=1/2;
F_end = 0;
L_sp=zeros(N/2+1,nScore,T);
F_sp=zeros(N/2+1,nScore,T,2);
L_team=zeros(N/2+1,nScore,T);
F_team=F_sp;
psi=1-lambda1-lambda2+2*lambda2*[0:1:N/2]'/N;
Psub_f_team=F_sp;
Psub_f_sp=F_sp;
Psub_l_team=L_sp;
Psub_l_sp=L_sp;
Pteam=F_sp;

L_team(:,:,end)=L_team_end;
L_sp(:,:,end)=L_sp_end;


%%

K = @(x,sigma) x.^sigma;
EC_Z = @(z,sigma) z*sigma/(1+sigma); %E[c|c<z] (truncated expectation)
fun_p = @(q,V1,V2,sigma) K(max(q*(V1-V2),0),sigma);
fun_exp_max = @(p,q,V1,V2,sigma) p*(q*V1+(1-q)*V2-EC_Z(q*(V1-V2),sigma))+(1-p)*V2;

for k=1:T-1
    t = T - k;
    for i=0:(N/2)
        aux_Teams = i;
        iaux_Teams = aux_Teams+1;
        aux_sp=N - 2*aux_Teams;
%         L_sp(:,nScore,t)=L_sp_end;
%         L_team(:,nScore,t)=L_team_end;
        for j=1:nScore-1
            aux_Score = j;
            
            p_team_f_sp_0=0;
            p_team_f_sp_1=0;
            F_sp_forms_team_0=0;
            F_sp_forms_team_1=0;
            F_sp_rival_sp_0=0;
            F_sp_rival_sp_1=0;
            F_sp_rival_team_0=0;
            F_sp_rival_team_1=0;
            F_sp_team_forms_0=0;
            F_sp_team_forms_1=0;
            L_sp_team_forms=0;
            L_team_team_forms=0;
            F_team_rival_team_0=0;
            F_team_rival_team_1=0;
            F_team_rival_sp_0=0;
            F_team_rival_sp_1=0;
            F_team_team_forms_0=0;
            F_team_team_forms_1=0;
            
             %%%%%%%%CCPs
            p_sub_f_sp_0 = fun_p(q_sp(aux_Score),L_sp(iaux_Teams,aux_Score+1,t+1),F_sp(iaux_Teams,aux_Score,t+1,1),sigmaS);
            p_sub_f_sp_1 = fun_p(q_sp(aux_Score),L_sp(iaux_Teams,aux_Score+1,t+1),F_sp(iaux_Teams,aux_Score,t+1,2),sigmaS);
            p_sub_f_team_0 = fun_p(q_team(aux_Score),L_team(iaux_Teams,aux_Score+1,t+1),F_team(iaux_Teams,aux_Score,t+1,1),sigmaS);
            p_sub_f_team_1 = fun_p(q_team(aux_Score),L_team(iaux_Teams,aux_Score+1,t+1),F_team(iaux_Teams,aux_Score,t+1,2),sigmaS);
            p_sub_l_sp_0 = fun_p(q_sp(aux_Score),L_sp(iaux_Teams,aux_Score+1,t+1),L_sp(iaux_Teams,aux_Score,t+1),sigmaS);
            p_sub_l_team_1 = fun_p(q_team(aux_Score),L_team(iaux_Teams,aux_Score+1,t+1),L_team(iaux_Teams,aux_Score,t+1),sigmaS);
            if aux_Teams<N/2
                p_team_f_sp_0 = fun_p(1,F_team(iaux_Teams+1,aux_Score,t+1,1),F_sp(iaux_Teams,aux_Score,t+1,1),sigmaT);
                p_team_f_sp_1 = fun_p(1,F_team(iaux_Teams+1,aux_Score,t+1,2),F_sp(iaux_Teams,aux_Score,t+1,2),sigmaT);
                
             %aux_F_sp_forms_team computes F_sp_forms team as a function of
            %F^{sp}_{s,n,l}, l, and p^team
                aux_F_sp_forms_team = @(p,l) fun_exp_max(p,1,F_team(iaux_Teams+1,aux_Score,t+1,1+l),F_sp(iaux_Teams,aux_Score,t+1,1+l),sigmaT);
                F_sp_forms_team_0 = aux_F_sp_forms_team(p_team_f_sp_0,0);
                F_sp_forms_team_1 = aux_F_sp_forms_team(p_team_f_sp_1,1);                
 
                %computing F^{sp,rival sp}
                aux_F_sp_rival_sp= @(pl,pf,F,l) ((1-l)/(aux_sp-1))*(pl*(q_sp(aux_Score)*F_sp(iaux_Teams,aux_Score+1,t+1,1+0)+(1-q_sp(aux_Score))*F_sp(iaux_Teams,aux_Score,t+1,1+0))+(1-pl)*F_sp(iaux_Teams,aux_Score,t+1,1+0))+((aux_sp-(1-l))/(aux_sp-1))*(pf*(q_sp(aux_Score)*F_sp(iaux_Teams,aux_Score+1,t+1,1+0)+(1-q_sp(aux_Score))*F)+(1-pf)*F);
                F_sp_rival_sp_0 = aux_F_sp_rival_sp(p_sub_l_sp_0,p_sub_f_sp_0,F_sp(iaux_Teams,aux_Score,t+1,1+0),0);
                F_sp_rival_sp_1 = aux_F_sp_rival_sp(p_sub_l_sp_0,p_sub_f_sp_1,F_sp(iaux_Teams,aux_Score,t+1,1+1),1);
                L_sp_team_forms = p_team_f_sp_0*L_sp(iaux_Teams+1,aux_Score,t+1) + (1-p_team_f_sp_0)*L_sp(iaux_Teams,aux_Score,t+1);
                L_team_team_forms = p_team_f_sp_1*L_team(iaux_Teams+1,aux_Score,t+1) + (1-p_team_f_sp_1)*L_team(iaux_Teams,aux_Score,t+1);
               
                aux_F_team_rival_sp= @(pl,pf,l) ((1-l)/(aux_sp))*(pl*(q_sp(aux_Score)*F_team(iaux_Teams,aux_Score+1,t+1,1+0)+(1-q_sp(aux_Score))*F_team(iaux_Teams,aux_Score,t+1,1+0))+(1-pl)*F_team(iaux_Teams,aux_Score,t+1,1+0))+((aux_sp-1+l)/(aux_sp))*(pf*(q_sp(aux_Score)*F_team(iaux_Teams,aux_Score+1,t+1,1+0)+(1-q_sp(aux_Score))*F_team(iaux_Teams,aux_Score,t+1,1+l))+(1-pf)*F_team(iaux_Teams,aux_Score,t+1,1+l));
                F_team_rival_sp_0 = aux_F_team_rival_sp(p_sub_l_sp_0,p_sub_f_sp_0,0);
                F_team_rival_sp_1 = aux_F_team_rival_sp(p_sub_l_sp_0,p_sub_f_sp_1,1);
                F_team_team_forms_0 = p_team_f_sp_0*F_team(iaux_Teams+1,aux_Score,t+1,1+0) + (1-p_team_f_sp_0)*F_team(iaux_Teams,aux_Score,t+1,1+0);
                F_team_team_forms_1 = p_team_f_sp_1*F_team(iaux_Teams+1,aux_Score,t+1,1+1) + (1-p_team_f_sp_1)*F_team(iaux_Teams,aux_Score,t+1,1+1);
                
            end
            
              
                
            %%%%%%%%%F_sp
            %aux_F_sp_own computes F_sp_own as a function of F^{sp}_{s,n,l} and
            %p^sub
            aux_F_sp_own = @(p,F) fun_exp_max(p,q_sp(aux_Score),L_sp(iaux_Teams,aux_Score+1,t+1),F,sigmaS);
            F_sp_own_0 = aux_F_sp_own(p_sub_f_sp_0,F_sp(iaux_Teams,aux_Score,t+1,1));
            F_sp_own_1 = aux_F_sp_own(p_sub_f_sp_1,F_sp(iaux_Teams,aux_Score,t+1,2));
            

            %computing F^{sp,rival team}
            if aux_Teams>0
                aux_F_sp_rival_team = @(pl,pf,F,l) (l/aux_Teams)*(pl*(q_team(aux_Score)*F_sp(iaux_Teams,aux_Score+1,t+1,1+1)+(1-q_team(aux_Score))*F_sp(iaux_Teams,aux_Score,t+1,1+1))+(1-pl)*F_sp(iaux_Teams,aux_Score,t+1,1+1))+((aux_Teams-l)/aux_Teams)*(pf*(q_team(aux_Score)*F_sp(iaux_Teams,aux_Score+1,t+1,1+1)+(1-q_team(aux_Score))*F)+(1-pf)*F);
                F_sp_rival_team_0 = aux_F_sp_rival_team(p_sub_l_team_1,p_sub_f_team_0,F_sp(iaux_Teams,aux_Score,t+1,1+0),0);
                F_sp_rival_team_1 = aux_F_sp_rival_team(p_sub_l_team_1,p_sub_f_team_1,F_sp(iaux_Teams,aux_Score,t+1,1+1),1);

            end
            
            %Computing F_sp_team forms
            if aux_sp>2
                aux_F_sp_team_forms = @(p,l) ((1-p)*((aux_sp-2+l)/(aux_sp-1))+((1-l)/(aux_sp-1)))*F_sp(iaux_Teams,aux_Score,t+1,1+l) + p*((aux_sp-2+l)/(aux_sp-1))*((1/(aux_sp-2+l))*F_team(iaux_Teams+1,aux_Score,t+1,1+l) + ((aux_sp-3+l)/(aux_sp-2+l))*F_sp(iaux_Teams+1,aux_Score,t+1,1+l));
                F_sp_team_forms_0 = aux_F_sp_team_forms(p_team_f_sp_0,0);
                F_sp_team_forms_1 = aux_F_sp_team_forms(p_team_f_sp_1,1);
            end
            
            F_sp(iaux_Teams,aux_Score,t,1) = (lambda1/N)*F_sp_own_0 + (lambda2/N)*F_sp_forms_team_0 + psi(iaux_Teams)*F_sp(iaux_Teams,aux_Score,t+1,1) + (lambda1*(2*aux_Teams)/N)*F_sp_rival_team_0 + lambda1*((aux_sp-1)/N)*F_sp_rival_sp_0 + lambda2*((aux_sp-1)/N)*F_sp_team_forms_0;
            F_sp(iaux_Teams,aux_Score,t,2) = (lambda1/N)*F_sp_own_1 + (lambda2/N)*F_sp_forms_team_1 + psi(iaux_Teams)*F_sp(iaux_Teams,aux_Score,t+1,2) + (lambda1*(2*aux_Teams)/N)*F_sp_rival_team_1 + lambda1*((aux_sp-1)/N)*F_sp_rival_sp_1 + lambda2*((aux_sp-1)/N)*F_sp_team_forms_1;
            
            %%%%%%%%%L_sp
            L_sp_own_0 = fun_exp_max(p_sub_l_sp_0,q_sp(aux_Score),L_sp(iaux_Teams,aux_Score+1,t+1),L_sp(iaux_Teams,aux_Score,t+1),sigmaS);
            L_sp_rival_team = p_sub_f_team_0*(q_team(aux_Score)*F_sp(iaux_Teams,aux_Score+1,t+1,1+1)+(1-q_team(aux_Score))*L_sp(iaux_Teams,aux_Score,t+1))+(1-p_sub_f_team_0)*L_sp(iaux_Teams,aux_Score,t+1);
            L_sp_rival_sp = p_sub_f_sp_0*(q_sp(aux_Score)*F_sp(iaux_Teams,aux_Score+1,t+1,1+0)+(1-q_sp(aux_Score))*L_sp(iaux_Teams,aux_Score,t+1))+(1-p_sub_f_sp_0)*L_sp(iaux_Teams,aux_Score,t+1);

            L_sp(iaux_Teams,aux_Score,t) = (lambda1/N)*L_sp_own_0 + (psi(iaux_Teams)+lambda2/N)*L_sp(iaux_Teams,aux_Score,t+1) + (lambda1*(2*aux_Teams)/N)*L_sp_rival_team + lambda1*((aux_sp-1)/N)*L_sp_rival_sp + lambda2*((aux_sp-1)/N)*L_sp_team_forms;
            
            
            %%%%%%%%%L_team
            L_team_own_1 = fun_exp_max(p_sub_l_team_1,q_team(aux_Score),L_team(iaux_Teams,aux_Score+1,t+1),L_team(iaux_Teams,aux_Score,t+1),sigmaS);
            L_team_rival_team = p_sub_f_team_1*(q_team(aux_Score)*F_team(iaux_Teams,aux_Score+1,t+1,1+1)+(1-q_team(aux_Score))*L_team(iaux_Teams,aux_Score,t+1))+(1-p_sub_f_team_1)*L_team(iaux_Teams,aux_Score,t+1);
            L_team_rival_sp = p_sub_f_sp_1*(q_sp(aux_Score)*F_team(iaux_Teams,aux_Score+1,t+1,1+0)+(1-q_sp(aux_Score))*L_team(iaux_Teams,aux_Score,t+1))+(1-p_sub_f_sp_1)*L_team(iaux_Teams,aux_Score,t+1);

            L_team(iaux_Teams,aux_Score,t) = ((2*lambda1)/N)*L_team_own_1 + psi(iaux_Teams)*L_team(iaux_Teams,aux_Score,t+1) + (lambda1*(2*(aux_Teams-1))/N)*L_team_rival_team + lambda1*(aux_sp/N)*L_team_rival_sp + lambda2*(aux_sp/N)*L_team_team_forms;
            
            
            %%%%%%%%%%F_team
            f_team_own_0 = fun_exp_max(p_sub_f_team_0,q_team(aux_Score),L_team(iaux_Teams,aux_Score+1,t+1),F_team(iaux_Teams,aux_Score,t+1,1),sigmaS);
            f_team_own_1 = fun_exp_max(p_sub_f_team_1,q_team(aux_Score),L_team(iaux_Teams,aux_Score+1,t+1),F_team(iaux_Teams,aux_Score,t+1,2),sigmaS);
            
            %computing F^{team,rival team}
            if aux_Teams>1
                aux_F_team_rival_team = @(pl,pf,l) (l/(aux_Teams-1))*(pl*(q_team(aux_Score)*F_team(iaux_Teams,aux_Score+1,t+1,1+1)+(1-q_team(aux_Score))*F_team(iaux_Teams,aux_Score,t+1,1+1))+(1-pl)*F_team(iaux_Teams,aux_Score,t+1,1+1))+((aux_Teams-l-1)/(aux_Teams-1))*(pf*(q_team(aux_Score)*F_team(iaux_Teams,aux_Score+1,t+1,1+1)+(1-q_team(aux_Score))*F_team(iaux_Teams,aux_Score,t+1,1+l))+(1-pf)*F_team(iaux_Teams,aux_Score,t+1,1+l));
                F_team_rival_team_0 = aux_F_team_rival_team(p_sub_l_team_1,p_sub_f_team_0,0);
                F_team_rival_team_1 = aux_F_team_rival_team(p_sub_l_team_1,p_sub_f_team_1,1);
            end
            
    
            F_team(iaux_Teams,aux_Score,t,1+0) = ((2*lambda1)/N)*f_team_own_0 + psi(iaux_Teams)*F_team(iaux_Teams,aux_Score,t+1,1+0) + (lambda1*(2*(aux_Teams-1))/N)*F_team_rival_team_0 + lambda1*((aux_sp)/N)*F_team_rival_sp_0 + lambda2*((aux_sp)/N)*F_team_team_forms_0;
            F_team(iaux_Teams,aux_Score,t,1+1) = ((2*lambda1)/N)*f_team_own_1 + psi(iaux_Teams)*F_team(iaux_Teams,aux_Score,t+1,1+1) + (lambda1*(2*(aux_Teams-1))/N)*F_team_rival_team_1 + lambda1*((aux_sp)/N)*F_team_rival_sp_1 + lambda2*((aux_sp)/N)*F_team_team_forms_1;
                       
            %%%% Storing CCPs
            Psub_f_team(iaux_Teams,aux_Score,t,1)=p_sub_f_team_0;
            Psub_f_team(iaux_Teams,aux_Score,t,2)=p_sub_f_team_1;
            Psub_f_sp(iaux_Teams,aux_Score,t,1)=p_sub_f_sp_0;
            Psub_f_sp(iaux_Teams,aux_Score,t,2)=p_sub_f_sp_1;
            Psub_l_team(iaux_Teams,aux_Score,t)=p_sub_l_team_1;
            Psub_l_sp(iaux_Teams,aux_Score,t)=p_sub_l_sp_0;
            Pteam(iaux_Teams,aux_Score,t,1)=p_team_f_sp_0;
            Pteam(iaux_Teams,aux_Score,t,2)=p_team_f_sp_1;
                   
        end
    end
[k T];    
end


%%
Out=zeros(nSim,4);
for k=1:nSim
    state=zeros(T+1,3);
    subs=zeros(T,1);
    mergers=zeros(T,1);
              %[s,n^teams,\ell] 
    state(1,:)=[1,0,0];
    for t=1:T
        rng(t+T*k)
        aux_sub=unifrnd(0,1)<=lambda1;
        aux_Teams = state(t,2);
        aux_sp = N - 2*aux_Teams;
        iaux_Teams=aux_Teams+1;
        aux_Score = state(t,1);
        aux_ell = state(t,3);
        if aux_sub==1
            aux_turn_sp=unifrnd(0,1)<=aux_sp/N;
            if aux_turn_sp==1
                if aux_ell==0
                    aux_turn_leader = unifrnd(0,1)<=1/aux_sp;
                    if aux_turn_leader==1
                        subs(t,1) = unifrnd(0,1)<=Psub_l_sp(iaux_Teams,aux_Score,t);
                    else
                        subs(t,1) = unifrnd(0,1)<=Psub_f_sp(iaux_Teams,aux_Score,t,aux_ell+1);                        
                    end     
                else
                    subs(t,1) = unifrnd(0,1)<=Psub_f_sp(iaux_Teams,aux_Score,t,aux_ell+1);
                end
                if subs(t,1)==1 && unifrnd(0,1)<=q_sp(aux_Score)
                    state(t+1,1) = min(aux_Score+1,nScore);
                    state(t+1,2) = state(t,2);
                    state(t+1,3) = 0;
                else
                    state(t+1,1) = state(t,1);
                    state(t+1,2) = state(t,2);
                    state(t+1,3) = state(t,3);                    
                end
            else
                if aux_ell==0
                    subs(t,1) = unifrnd(0,1)<=Psub_f_team(iaux_Teams,aux_Score,t,aux_ell+1);
                else
                    aux_turn_leader = unifrnd(0,1)<=1/aux_Teams;                    
                    if aux_turn_leader==1
                        subs(t,1) = unifrnd(0,1)<=Psub_l_team(iaux_Teams,aux_Score,t);
                    else
                        subs(t,1) = unifrnd(0,1)<=Psub_f_team(iaux_Teams,aux_Score,t,aux_ell+1);                        
                    end                       
                end
                if subs(t,1)==1 && unifrnd(0,1)<=q_team(aux_Score)
                    state(t+1,1) = min(aux_Score+1,nScore);
                    state(t+1,2) = state(t,2);
                    state(t+1,3) = 1;
                else
                    state(t+1,1) = state(t,1);
                    state(t+1,2) = state(t,2);
                    state(t+1,3) = state(t,3);                    
                end
            end
        else
            if aux_Teams==N/2
                state(t+1,1) = state(t,1);
                state(t+1,2) = state(t,2);
                state(t+1,3) = state(t,3);
            else
                aux_turn_sp=unifrnd(0,1)<=(aux_sp-(1-aux_ell))/N;
                if aux_turn_sp==1
                    mergers(t,1) = unifrnd(0,1)<=Pteam(iaux_Teams,aux_Score,t);
                    state(t+1,1) = state(t,1);
                    state(t+1,2) = state(t,2)+mergers(t,1);
                    state(t+1,3) = state(t,3);
                else
                    state(t+1,1) = state(t,1);
                    state(t+1,2) = state(t,2);
                    state(t+1,3) = state(t,3);
                end
            end
            
        end
    end
    Out(k,:)=[sum(subs),Svalues(state(T,1)),sum(mergers),(sum(mergers)-mergers_data).^2];
end
Out_aux=mean(Out,1);
subs=Out_aux(1);
maxs=Out_aux(2);
mergers=Out_aux(3);
Omega=Out_aux(4);
end
